# R script for replicating results in appendix section "Robustness Checks"
# Zhai & Garside 2022
# original device: MacPro 13, R 4.1.2
# recommended working directory: Desktop

setwd("~/Desktop") # default wkd
rm(list = ls())

# pkgs --------------------------------------------------------------------

if (!require("pacman")) install.packages("pacman")
pacman::p_load(tidyverse, magrittr, fixest, knitr, texreg)

# helpers -----------------------------------------------------------------

source("R_00_helpers.R")

# table 18 ----------------------------------------------------------------

## data ----
df.main <- read.csv("data_vote_main.csv", row.names = 1) %>% 
  group_by(land_name) %>% 
  mutate(flooded_land = max(flooded, na.rm = TRUE)) %>% 
  ungroup() %>% 
  filter(flooded_land == 1) 

## models ----
fm.h1 <- "v_green_pct~factor(flooded)*factor(date)+pop_per_sqkm+net_out_per_thousand+old_share+income_mean+unemployed_rate+land_set_pct+land_agri_pct+factor(incumbent)+distance*factor(date)"
fm.h1.xx <- gsub("v_green_pct\\~factor\\(flooded\\)\\*factor\\(date\\)\\+","",fm.h1)
fm.h1.xx.int <- unlist(strsplit(fm.h1.xx, "\\+")) %>% paste0("*factor(flooded)") %>% paste(collapse = "+")
fm.h1.int <- paste0("v_green_pct~factor(flooded)*factor(date)+",fm.h1.xx.int)
fm2.h1.int <- gsub("flooded", "flooded2", fm.h1.int)
tab18 <- map(
  list(fm.h1.int,fm2.h1.int),
  ~fit_twfe(df = df.main, fml = .x, fe = "gen_id", cl = "gen_id"))

## table ----
screenreg(tab18, digits = 3)

# table 19 & figure 19 ----------------------------------------------------

## data ----
df.allparty.17 <- readRDS("data_vote_all_ts.RDS") %>% 
  filter(date==max(date))
df.allparty.21 <- readRDS("data_vote_all_main.RDS")
id.pct <- paste0("v_", c("cducsu", "spd", "afd", "fdp", "dielinke", "green"), "_pct")
df.allparty <- 
  rbind.data.frame(
    df.allparty.17 %>% select(all_of(id.pct), date, gen_id), 
    df.allparty.21 %>% select(all_of(id.pct), date, gen_id)
  ) %>% 
  left_join(df.main %>% select(-v_green_pct) %>% mutate(date=as.Date(date))) %>% 
  filter(land_name %in% df.main$land_name)

## models ----
fm.rhs.sur <- "~factor(flooded)*factor(date)+pop_per_sqkm+net_out_per_thousand+old_share+income_mean+unemployed_rate+land_set_pct+land_agri_pct+factor(incumbent)+distance*factor(date)"
fm.rhs.sur.x <- gsub("flooded", "flooded2", fm.rhs.sur)
fm.lhs.sur <- id.pct
fm.sur <- vector("list", length = length(fm.lhs.sur))
names(fm.sur) <- gsub("\\_","", fm.lhs.sur)
for (i in seq_along(fm.lhs.sur)) {
  fm.sur[[i]] <- as.formula(paste0(fm.lhs.sur[i], fm.rhs.sur))
}
fm.sur.x <- vector("list", length = length(fm.lhs.sur))
names(fm.sur.x) <- gsub("\\_","", fm.lhs.sur)
for (i in seq_along(fm.lhs.sur)) {
  fm.sur.x[[i]] <- as.formula(paste0(fm.lhs.sur[i], fm.rhs.sur.x))
}

## table 19 ----
pacman::p_load(systemfit) # avoid ^namespace conflict (MASS::select vs dplyr::select)
tab19 <- 
  list(
    systemfit(formula = fm.sur, method = "SUR", data = df.allparty),
    systemfit(formula = fm.sur.x, method = "SUR", data = df.allparty)
  )
pacman::p_unload(systemfit) 
select <- dplyr::select # force namespace clear (select = dplyr::select)

screenreg(tab19, digits = 3, 
          custom.coef.map = list(
            "vcducsupct: factor(flooded)1:factor(date)2021-09-26"="CDU/CSU",
            "vspdpct: factor(flooded)1:factor(date)2021-09-26"="SPD",
            "vafdpct: factor(flooded)1:factor(date)2021-09-26"="AfD",
            "vfdppct: factor(flooded)1:factor(date)2021-09-26"="FDP",
            "vdielinkepct: factor(flooded)1:factor(date)2021-09-26"="Die Linke",
            "vgreenpct: factor(flooded)1:factor(date)2021-09-26"="Greens",
            "vcducsupct: factor(flooded2)1:factor(date)2021-09-26"="CDU/CSU",
            "vspdpct: factor(flooded2)1:factor(date)2021-09-26"="SPD",
            "vafdpct: factor(flooded2)1:factor(date)2021-09-26"="AfD",
            "vfdppct: factor(flooded2)1:factor(date)2021-09-26"="FDP",
            "vdielinkepct: factor(flooded2)1:factor(date)2021-09-26"="Die Linke",
            "vgreenpct: factor(flooded2)1:factor(date)2021-09-26"="Greens"))

## figure 19 ----
ls.fig19 <- vector("list", length = length(tab19))
for (i in seq_along(tab19)) {
  ls.fig19[[i]] <- 
    map(
      list(tab19[[i]]$coefficients, tab19[[i]]$coefCov %>% diag() %>% sqrt()),
      ~{.x[grepl("factor\\((flooded|flooded2)\\)1\\:factor\\(date\\)2021-09-26",names(.x))]}) %>%
    map(~purrr::set_names(.x, fm.lhs.sur)) %>%
    `names<-`(c("mu","se")) %>%
    bind_rows(.id = "measure") %>%
    pivot_longer(cols = -measure, names_to = "party", values_to = "est") %>%
    pivot_wider(names_from = measure, values_from = est) %>%
    mutate(lwr = mu-1.96*se, upr = mu+1.96*se) %>%
    mutate(party = case_when(
      party == "v_cducsu_pct" ~"CDU/CSU",
      party == "v_spd_pct" ~ "SPD",
      party == "v_afd_pct" ~ "AfD",
      party == "v_fdp_pct" ~ "FDP",
      party == "v_dielinke_pct" ~ "Die Linke",
      party == "v_green_pct" ~ "Greens"))
}  
names(ls.fig19) <- c("Primary","Secondary")
df.fig19 <- ls.fig19 %>% bind_rows(.id = "measure")

ggplot(df.fig19) +
  geom_pointrange(aes(x=party, y=mu, ymin=lwr, ymax=upr, color = party, shape = party)) +
  facet_grid(measure~.) +
  geom_hline(yintercept = 0, color = "darkred", lty = "dotted") +
  coord_flip() +
  ggthemes::scale_color_stata() +
  labs(x = NULL, y = NULL, color = "Party", shape = "Party") +
  scale_y_continuous(breaks = scales::pretty_breaks(n=5)) +
  theme_minimal(base_size = 14) 
ggsave("appendix_figure19.pdf", width = 9, height = 5)

# figure 20 ---------------------------------------------------------------

## data ----
IDs.state <- unique(df.main$land_name)
df.allparty.state <- make_state_dfs(df.allparty)

## models ----
fm.sur.state <- drop_incumbent(fml = fm.sur)
fm.sur.x.state <- drop_incumbent(fml = fm.sur.x)

ls.fig20 <- replicate(2, vector("list", length = length(IDs.state)), simplify = FALSE)
names(ls.fig20[[1]]) <- names(ls.fig20[[2]]) <- IDs.state
for (j in seq_along(IDs.state)) {
  ls.fig20[[1]][[j]] <- sfit_sur(df = df.allparty.state[[j]], fml = fm.sur.state)
  ls.fig20[[2]][[j]] <- sfit_sur(df = df.allparty.state[[j]], fml = fm.sur.x.state)
}

## plot ----
plot_sur_est(ls.fig20, lhs = fm.lhs.sur)
ggsave("~/Desktop/appendix_figure20.pdf", width = 18, height = 5)

# table 20 ----------------------------------------------------------------

df.main0 <- read.csv("data_vote_main.csv", row.names = 1)
xx.h1 <- "+pop_per_sqkm+net_out_per_thousand+old_share+income_mean+unemployed_rate+land_set_pct+land_agri_pct+factor(incumbent)+distance*factor(date)"
fm.h1 <- "v_green_pct~factor(flooded)*factor(date)"
fm.h1.x <- paste0(fm.h1, xx.h1)
fm2.h1 <- gsub("flooded", "flooded2", fm.h1)
fm2.h1.x <- paste0(fm2.h1, xx.h1)

tab20 <- map(
  list(fm.h1, fm.h1.x, fm2.h1, fm2.h1.x),
  ~fit_twfe(df = df.main0, fm = .x, fe = "gen_id", cl = "gen_id"))
screenreg(tab20, digits = 3)

# table 21 ----------------------------------------------------------------

df.main0.h1a <- df.main0 %>% 
  mutate(nsevere = as.integer(flooded==1 & severe==0),
         nsevere2 = as.integer(flooded2==1 & severe2==0))
fm.h1a <- "v_green_pct~factor(severe)*factor(date)+factor(nsevere)*factor(date)"
fm.h1a.x <- paste0(fm.h1a, xx.h1)
fm2.h1a <- gsub("severe", "severe2", fm.h1a)
fm2.h1a.x <- paste0(fm2.h1a, xx.h1)

tab21 <- map(
  list(fm.h1a, fm.h1a.x, fm2.h1a, fm2.h1a.x),
  ~fit_twfe(df = df.main0.h1a, fm = .x, fe = "gen_id", cl = "gen_id"))
screenreg(tab21, digits = 3)

# table 22 ----------------------------------------------------------------

## data----
df.main0 %<>% group_by(land_name) %>% 
  mutate(flooded_land = max(flooded, na.rm = TRUE)) %>% 
  ungroup()
df.spatmat <- readRDS("data_spatmat.RDS") 
df.spatmat.knn <- 
  df.spatmat$knn %>% 
  as.data.frame() %>% 
  mutate(gen_id = as.integer(colnames(.)))
var.h1b <- 
  df.main0 %>% 
  filter(date == max(date)) %>% 
  select(gen_id, land_name, flooded, flooded2) %>% 
  left_join(df.spatmat.knn, by = "gen_id")
d.h1b <- as.matrix(var.h1b[,names(var.h1b) %in% df.spatmat.knn$gen_id]) %*% var.h1b$flooded 
d2.h1b <- as.matrix(var.h1b[,names(var.h1b) %in% df.spatmat.knn$gen_id]) %*% var.h1b$flooded2
df.h1b <- 
  data.frame(
    gen_id = rep(var.h1b$gen_id, 2),
    flooded_w = rep(d.h1b, 2),
    flooded2_w = rep(d2.h1b, 2),
    date = c(rep("2017-09-24", length(d2.h1b)), rep("2021-09-26", length(d2.h1b)))
  ) %>% 
  left_join(df.main0)

## models ----
fm.h1b <- "v_green_pct~factor(flooded2)*factor(date)+flooded2_w*factor(date)"
fm.h1b.x <- paste0(fm.h1b, xx.h1)

## table ----
tab22 <- map(
  list(fm.h1b, fm.h1b.x),
  ~fit_twfe(df = df.h1b, fm = .x, fe = "gen_id", cl = "gen_id"))
screenreg(tab22, digits = 3)

# tables 23 & 24 ----------------------------------------------------------

## data ----
df.pcpt0 <- read.csv("data_gles_wknr.csv", row.names = 1)
df.pcpt <- 
  df.pcpt0 %>% 
  select(wk_id, year, climate_position_pre = climate_position, climate_salience_pre = climate_salience) %>% 
  filter(year==min(year)) %>% 
  select(-year) %>% 
  left_join(df.pcpt0, ., by = "wk_id") 

## models ----
xx.pcpt <- "+pop_per_sqkm+net_out_per_thousand+old_share+income_mean+unemployed_rate+land_set_pct+land_agri_pct+factor(incumbent)+distance*factor(year)"
fm.pcpt.a <- "v_green_pct~factor(flooded)*factor(year)*climate_position_pre"
fm.pcpt.b <- "v_green_pct~factor(flooded)*factor(year)*climate_salience_pre"
fm.pcpt.x.a <- paste0(fm.pcpt.a, xx.pcpt)
fm.pcpt.x.b <- paste0(fm.pcpt.b, xx.pcpt)
fm2.pcpt.a <- gsub("flooded", "flooded2", fm.pcpt.a)
fm2.pcpt.b <- gsub("flooded", "flooded2", fm.pcpt.b)
fm2.pcpt.x.a <- paste0(fm2.pcpt.a, xx.pcpt)
fm2.pcpt.x.b <- paste0(fm2.pcpt.b, xx.pcpt)

## table 23 ----
tab23 <- map(
  list(fm.pcpt.a, fm.pcpt.x.a, fm2.pcpt.a, fm2.pcpt.x.a), 
  ~fit_twfe(df = df.pcpt, fm = .x, fe = "wk_id", cl = "wk_id"))
screenreg(tab23, digits = 3)

## table 24 ----
tab24 <- map(
  list(fm.pcpt.b, fm.pcpt.x.b, fm2.pcpt.b, fm2.pcpt.x.b),
  ~fit_twfe(df = df.pcpt, fm = .x, fe = "wk_id", cl = "wk_id"))
screenreg(tab24, digits = 3)

# table 25 ----------------------------------------------------------------

xx.h2 <- "+pop_per_sqkm+net_out_per_thousand+old_share+income_mean+unemployed_rate+land_set_pct+land_agri_pct+factor(incumbent)+distance*factor(year)"
fm.h2.a <- "climate_position~factor(flooded)*factor(year)"
fm.h2.b <- "climate_salience~factor(flooded)*factor(year)"
fm.h2.x.a <- paste0(fm.h2.a, xx.h2)
fm.h2.x.b <- paste0(fm.h2.b, xx.h2)
fm2.h2.a <- gsub("flooded", "flooded2", fm.h2.a)
fm2.h2.b <- gsub("flooded", "flooded2", fm.h2.b)
fm2.h2.x.a <- paste0(fm2.h2.a, xx.h2)
fm2.h2.x.b <- paste0(fm2.h2.b, xx.h2)

tab25.a <- map(
  list(fm.h2.a, fm.h2.x.a, fm2.h2.a, fm2.h2.x.a),
  ~fit_twfe(df = df.pcpt0, fm = .x, fe = "wk_id", cl = "wk_id"))
tab25.b <- map(
  list(fm.h2.b, fm.h2.x.b, fm2.h2.b, fm2.h2.x.b),
  ~fit_twfe(df = df.pcpt0, fm = .x, fe = "wk_id", cl = "wk_id"))
screenreg(c(tab25.a, tab25.b), digits = 3, custom.header = list("Climate Position"=1:4, "Issue Salience"=5:8))

# table 26 ----------------------------------------------------------------

df.h3a <- readRDS("data_turnout_main.RDS") %>% 
  left_join(df.main0 %>% filter(date == max(date)) %>% select(gen_id, starts_with("flood"), v_green_pct))

tab26 <- 
  rbind.data.frame(
    df.h3a %>% nest(-flooded) %>% 
      mutate(corr = map(data, ~cor.test(.$turnout, .x$v_green_pct)) %>% map(., broom::tidy)) %>% 
      unnest(corr) %>% select(-data) %>% 
      add_column(measure = "primary"),
    df.h3a %>% nest(-flooded2) %>% 
      mutate(corr = map(data, ~cor.test(.$turnout, .x$v_green_pct)) %>% map(., broom::tidy)) %>% 
      unnest(corr) %>% select(-data) %>% 
      rename(flooded = flooded2) %>% 
      add_column(measure = "secondary"))
tab26[,c("flooded","estimate","statistic","p.value","measure")] %>% knitr::kable()

# table 27 ----------------------------------------------------------------

df.h3b <- read.csv("~/Dropbox/Projects/german-flood/7_submission/rap/04_article/replication_files/data_gles_individual.csv", row.names = 1) %>% 
  mutate(v_change = v_2021-v_2017)
xx.h3b <- "+factor(gender)+age+factor(married)+factor(immigrant)+educ+factor(empl)+hhinc+urban+relig+lr"
fm.h3b <- "v_change~factor(flooded)"
fm.h3b.x <- paste0(fm.h3b, xx.h3b)
fm2.h3b <- "v_change~factor(flooded2)"
fm2.h3b.x <- paste0(fm2.h3b, xx.h3b)

tab27 <- map(
  list(fm.h3b, fm.h3b.x, fm2.h3b, fm2.h3b.x),
  ~fit_twfe2(df = df.h3b, fml = .x, cl = "wk_id", wt = df.h3b$w_ow))
screenreg(tab27, digits = 3)

# table 28 ----------------------------------------------------------------

fm.h1.xx <- gsub("v_green_pct\\~factor\\(flooded\\)\\*factor\\(date\\)\\+","",fm.h1.x)
fm.h1.xx.int <- unlist(strsplit(fm.h1.xx, "\\+")) %>% paste0("*factor(flooded)") %>% paste(collapse = "+")
fm.h1.int <- paste0("v_green_pct~factor(flooded)*factor(date)+",fm.h1.xx.int)
fm2.h1.int <- gsub("flooded", "flooded2", fm.h1.int)

tab28 <- map(list(fm.h1.int,fm2.h1.int), 
             ~fit_twfe(df = df.main0, fml = .x, fe = "gen_id", cl = "gen_id"))
screenreg(tab28, digits = 3)

# table 29 ----------------------------------------------------------------

## data ----
df.med.turnout <- rbind.data.frame(
  readRDS("data_turnout_main.RDS") %>% 
    select(gen_id, date, turnout, date),
  readRDS("data_turnout_ts.RDS") %>% 
    filter(date==max(date)) %>% 
    mutate(turnout = turnout/100) %>% 
    select(gen_id, date, turnout)
)
df.med.turnout.fd <- 
  df.med.turnout %>% 
  mutate(date=as.Date(date), gen_id = as.integer(gen_id)) %>% 
  group_by(gen_id) %>% 
  summarise(turnout = turnout - lag(turnout)) %>% 
  ungroup() %>% 
  na.omit()
df.med.fd <- 
  df.main0 %>% 
  mutate(date=as.Date(date), gen_id = as.integer(gen_id)) %>% 
  group_by(gen_id) %>% 
  summarise(across(.cols = c(v_green_pct, pop_per_sqkm:incumbent),
                   .fns = ~{.x - lag(.x)})) %>% 
  ungroup() %>% 
  na.omit()
df.med.turnouts <- 
  left_join(df.med.fd, df.med.turnout.fd, by = "gen_id") %>% 
  left_join(
    df.main0 %>% dplyr::select(gen_id, flooded, flooded2) %>% distinct(),
    by = "gen_id"
  ) %>% 
  mutate(
    flooded=factor(flooded), 
    flooded2=factor(flooded2))

## models ----
pacman::p_load(mediation) # avoid ^namespace conflict (MASS::select vs dplyr::select)
set.seed(1)

tab29.med <- lm(turnout~flooded+pop_per_sqkm+net_out_per_thousand+old_share+income_mean+unemployed_rate+land_set_pct+land_agri_pct+factor(incumbent), data = df.med.turnouts)
tab29.out <- lm(v_green_pct~turnout+flooded+pop_per_sqkm+net_out_per_thousand+old_share+income_mean+unemployed_rate+land_set_pct+land_agri_pct+factor(incumbent), data = df.med.turnouts)
tab29.1 <- mediate(tab29.med, tab29.out, treat = "flooded", mediator = "turnout", robustSE = TRUE)
tab29.med2 <- lm(turnout~flooded2+pop_per_sqkm+net_out_per_thousand+old_share+income_mean+unemployed_rate+land_set_pct+land_agri_pct+factor(incumbent), data = df.med.turnouts)
tab29.out2 <- lm(v_green_pct~turnout+flooded2+pop_per_sqkm+net_out_per_thousand+old_share+income_mean+unemployed_rate+land_set_pct+land_agri_pct+factor(incumbent), data = df.med.turnouts)
tab29.2 <- mediate(tab29.med2, tab29.out2, treat = "flooded2", mediator = "turnout", robustSE = TRUE)
tab29 <- list(tab29.1,tab29.2)

pacman::p_unload(mediation)
select <- dplyr::select # force namespace clear (select = dplyr::select)

## summaries ----
summary(tab29.1)
summary(tab29.2)

# table 30 ----------------------------------------------------------------

## data ----
df.allparty0 <- 
  rbind.data.frame(
    df.allparty.17 %>% select(all_of(id.pct), date, gen_id), 
    df.allparty.21 %>% select(all_of(id.pct), date, gen_id)
  ) %>% 
  left_join(df.main %>% select(-v_green_pct) %>% mutate(date=as.Date(date)))

## models ----
pacman::p_load(systemfit) # avoid ^namespace conflict (MASS::select vs dplyr::select)
tab30 <- 
  list(
    systemfit(formula = fm.sur, method = "SUR", data = df.allparty0),
    systemfit(formula = fm.sur.x, method = "SUR", data = df.allparty0)
  )
pacman::p_unload(systemfit) 
select <- dplyr::select # force namespace clear (select = dplyr::select)

## table ----
screenreg(tab30, digits = 3, 
          custom.coef.map = list(
            "vcducsupct: factor(flooded)1:factor(date)2021-09-26"="CDU/CSU",
            "vspdpct: factor(flooded)1:factor(date)2021-09-26"="SPD",
            "vafdpct: factor(flooded)1:factor(date)2021-09-26"="AfD",
            "vfdppct: factor(flooded)1:factor(date)2021-09-26"="FDP",
            "vdielinkepct: factor(flooded)1:factor(date)2021-09-26"="Die Linke",
            "vgreenpct: factor(flooded)1:factor(date)2021-09-26"="Greens",
            "vcducsupct: factor(flooded2)1:factor(date)2021-09-26"="CDU/CSU",
            "vspdpct: factor(flooded2)1:factor(date)2021-09-26"="SPD",
            "vafdpct: factor(flooded2)1:factor(date)2021-09-26"="AfD",
            "vfdppct: factor(flooded2)1:factor(date)2021-09-26"="FDP",
            "vdielinkepct: factor(flooded2)1:factor(date)2021-09-26"="Die Linke",
            "vgreenpct: factor(flooded2)1:factor(date)2021-09-26"="Greens"))

# table 31 ----------------------------------------------------------------

df.gm <- readRDS("data_genmatch_main.RDS") # genetically matched datasets (main)

tab31 <- map2(df.gm[1:2], list(fm.h1.x, fm2.h1.x),~fit_twfe(.x, fml = .y, fe = "gen_id", cl = "gen_id"))
screenreg(tab31, digits = 3)

# table 32 ----------------------------------------------------------------

df.h1a.gm <- map(
  df.gm[grepl("severe", names(df.gm))],
  ~mutate(.x, nsevere = as.integer(flooded==1 & severe==0), nsevere2 = as.integer(flooded2==1 & severe2==0)))

tab32 <- fit_twfe(df.h1a.gm$severe, fm.h1a.x, fe = "gen_id", cl = "gen_id")
screenreg(tab32)

# table 33 ----------------------------------------------------------------

d2.h1b.gm <- as.matrix(df.spatmat.knn[as.character(df.gm$flooded2$gen_id),as.character(df.gm$flooded2$gen_id)]) %*% df.gm$flooded2$flooded2
df2.h1b.gm <- cbind.data.frame(df.gm$flooded2, "flooded2_w" = d2.h1b.gm)

tab33 <- fit_twfe(df2.h1b.gm, fm.h1b.x, fe="gen_id", cl="gen_id")
screenreg(tab33, digits = 3)

# table 34 ----------------------------------------------------------------

## data ----
df.gm.wk <- readRDS("data_genmatch_wknr.RDS") # genmatched datasets (Wahlkreis-level)

## models ----
fm.h2.x.a.gm <- "climate_position~factor(flooded)*factor(year)+pop_per_sqkm+net_out_per_thousand+old_share+income_mean+unemployed_rate+land_set_pct+land_agri_pct+factor(incumbent)+distance*factor(year)"
fm.h2.x.b.gm <- gsub("climate_position", "climate_salience", fm.h2.x.a.gm)
fm2.h2.x.a.gm <- 
  gsub("flooded", "flooded2", fm.h2.x.a.gm) %>% 
  gsub("\\+factor\\(incumbent\\)", "", .) # incumbent = 0 \forall wk
fm2.h2.x.b.gm <- gsub("climate_position", "climate_salience", fm2.h2.x.a.gm)

## table ----

tab34.a <- map2(
  list(fm.h2.x.a.gm, fm2.h2.x.a.gm), df.gm.wk,
  ~fit_twfe(df = .y, fml = .x, fe = "wk_id", cl = "wk_id"))
tab34.b <- map2(
  list(fm.h2.x.b.gm, fm2.h2.x.b.gm), df.gm.wk,
  ~fit_twfe(df = .y, fml = .x, fe = "wk_id", cl = "wk_id"))
screenreg(c(tab34.a, tab34.b), digits = 3, custom.header = list("Climate Position"=1:2, "Issue Salience"=3:4))

# table 35 ----------------------------------------------------------------

## data ----
df.h3a.gm <- 
  df.h3a %>% 
  filter(gen_id %in% df.gm$flooded$gen_id) %>% 
  left_join(df.gm$flooded %>% filter(date == max(date)) %>% select(gen_id, flooded, v_green_pct))
df2.h3a.gm <- 
  df.h3a %>% 
  filter(gen_id %in% df.gm$flooded2$gen_id) %>% 
  left_join(df.gm$flooded2 %>% filter(date == max(date)) %>% select(gen_id, flooded2, v_green_pct))

## table
tab35 <- rbind.data.frame(
  test_h3a(df.h3a.gm, D = flooded, measure = "primary"),
  test_h3a(df2.h3a.gm, D = flooded2, measure = "secondary")
)
tab35[,c("flooded","estimate","statistic","p.value","measure")] %>% knitr::kable()

# table 36 ----------------------------------------------------------------

## data ----
df.h3b.gm <- df.h3b %>% 
  filter(wk_id %in% df.gm.wk$flooded$wk_id)
df2.h3b.gm <- df.h3b %>% 
  filter(wk_id %in% df.gm.wk$flooded2$wk_id)

## table ----
tab36 <- map2(
  list(fm.h3b.x, fm2.h3b.x), list(df.h3b.gm, df2.h3b.gm),
  ~fit_twfe2(df = .y, fml = .x, cl = "wk_id", wt = .y[["w_ow"]]))
screenreg(tab36, digits = 3)

# table 37 ----------------------------------------------------------------

df.risk <- readRDS("data_flood_risk.RDS") # flood risk datasets
df.risk.main <- df.main0 %>% 
  select(-gen_name) %>% 
  left_join(df.risk[[1]], ., by ="gen_id")

tab37 <- map(
  list(fm.h1.x, fm2.h1.x),
  ~fit_twfe(df = df.risk.main, fml = .x, fe = "gen_id", cl = "gen_id"))
screenreg(tab37, digits = 3)

# table 38 ----------------------------------------------------------------

df.h1a.risk <- df.risk.main %>% 
  mutate(nsevere = as.integer(flooded==1 & severe==0),
         nsevere2 = as.integer(flooded2==1 & severe2==0))

tab38 <- map(
  list(fm.h1a.x, fm2.h1a.x),
  ~fit_twfe(df = df.h1a.risk, fml = .x, fe = "gen_id", cl = "gen_id")) 
screenreg(tab38, digits = 3)

# table 39 ----------------------------------------------------------------

d2.h1b.risk <- as.matrix(df.spatmat.knn[as.character(df.risk[[1]]$gen_id),as.character(df.risk[[1]]$gen_id)]) %*% df.risk.main$flooded2[df.risk.main$date==max(df.risk.main$date)]
df.h1b.risk <- cbind.data.frame(df.risk.main, "flooded2_w" = rep(d2.h1b.risk, each=2))

tab39 <- fit_twfe(df.h1b.risk, fm.h1b.x, fe="gen_id", cl="gen_id")
screenreg(tab39, digits = 3)

# table 40 ----------------------------------------------------------------

df.h2.risk <- df.pcpt0 %>% filter(wk_id %in% df.risk[[2]]$wk_id)

tab40.a <- map(
  list(fm.h2.x.a, fm2.h2.x.a),
  ~fit_twfe(df = df.h2.risk, fm = .x, fe = "wk_id", cl = "wk_id"))
tab40.b <- map(
  list(fm.h2.x.b, fm2.h2.x.b),
  ~fit_twfe(df = df.h2.risk, fm = .x, fe = "wk_id", cl = "wk_id"))
screenreg(c(tab40.a, tab40.b), digits = 3, custom.header = list("Climate Position"=1:2, "Issue Salience"=3:4))

# table 41 ----------------------------------------------------------------

df.h3a.risk <- df.h3a %>% 
  filter(gen_id %in% df.risk[[1]]$gen_id) %>% 
  left_join(df.risk.main %>% filter(date == max(date)) %>% select(gen_id,starts_with("flooded"), v_green_pct))

tab41 <- rbind.data.frame(
  test_h3a(df.h3a.risk, measure = "primary"),
  test_h3a(df.h3a.risk, D = flooded2, measure = "secondary"))
tab41[,c("flooded","estimate","statistic","p.value","measure")] %>% knitr::kable()

# table 42 ----------------------------------------------------------------

df.h3b.risk <- df.h3b %>% filter(wk_id %in% df.risk[[2]]$wk_id)

tab42 <- map(
  list(fm.h3b.x, fm2.h3b.x), 
  ~fit_twfe2(df = df.h3b.risk, fml = .x, cl = "wk_id", wt = df.h3b.risk[["w_ow"]]))
screenreg(tab42, digits = 3)

# table 43 ----------------------------------------------------------------

df.fv <- 
  readRDS("data_vote_1st.RDS") %>% 
  select(-gen_name) %>% 
  left_join(df.main0 %>% select(-starts_with("v_")) %>% mutate(date=as.Date(date))) %>% 
  select(gen_id, gen_name, date, everything()) 

tab43 <- map(
  list(fm.h1.x, fm2.h1.x),
  ~fit_twfe(df = df.fv, fm = .x, fe = "gen_id", cl = "gen_id"))
screenreg(tab43, digits = 3)

# table 44 ----------------------------------------------------------------

df.h1a.fv <- df.fv %>% 
  mutate(nsevere = as.integer(flooded==1 & severe==0),
         nsevere2 = as.integer(flooded2==1 & severe2==0))

tab44 <- map(
  list(fm.h1a.x, fm2.h1a.x), 
  ~fit_twfe(df = df.h1a.fv, fm = .x, fe = "gen_id", cl = "gen_id"))
screenreg(tab44, digits = 3)

# table 45 ----------------------------------------------------------------

d2.h1b.fv <- as.matrix(df.spatmat.knn[as.character(df.fv$gen_id[df.fv$date==max(df.fv$date)]),as.character(df.fv$gen_id[df.fv$date==max(df.fv$date)])]) %*% 
  df.fv$flooded2[df.fv$date==max(df.fv$date)]
df.h1b.fv <- cbind.data.frame(df.fv, "flooded2_w" = rep(d2.h1b.fv, 2))

tab45 <- fit_twfe(df = df.h1b.fv, fm = fm.h1b.x, fe = "gen_id", cl = "gen_id")
screenreg(tab45, digits = 3)

# table 46 ----------------------------------------------------------------

df.h3a.fv <- readRDS("data_turnout_1st.RDS") %>% 
  left_join(select(filter(df.fv, date==max(date)), gen_id, flooded, flooded2, v_green_pct))

tab46 <- rbind.data.frame(
    test_h3a(df.h3a.fv, measure = "primary"),
    test_h3a(df.h3a.fv, D = flooded2, measure = "secondary"))
tab46[,c("flooded","estimate","statistic","p.value","measure")] %>% knitr::kable()

# table 47 ----------------------------------------------------------------

df.h3b.fv <- read.csv("data_gles_1st.csv") %>% 
  mutate(v_change = v_2021-v_2017) 

tab47 <- map(
  list(fm.h3b.x, fm2.h3b.x), 
  ~fit_twfe2(df = df.h3b.fv, fml = .x, cl = "wk_id", wt = df.h3b.fv$w_ow))
screenreg(tab47, digits = 3)

# table 48 ----------------------------------------------------------------

## data ----
df.spatmat.queen <- as.data.frame(df.spatmat$queen)
d2.h1b.queen <- as.matrix(df.spatmat.queen[as.character(df.main0$gen_id),as.character(df.main0$gen_id)]) %*% df.main0$flooded2
df.queen <- cbind.data.frame(df.main0, "flooded2_w" = d2.h1b.queen) %>% 
  filter(gen_id %in% df.main$gen_id)

## table ----
tab48 <- fit_twfe(df = df.queen, fml = fm.h1b.x, fe = "gen_id", cl = "gen_id")
screenreg(tab48, digits = 3)

# cleanup -----------------------------------------------------------------

pacman::p_unload("all")
rm(list = ls())
